function [sol,mom,cal] = calibrate_production_function(version_pf,version_skill_dist,knitro_path)
%% paths
addpath(knitro_path)
addpath('../Densities')
addpath('../UtilityFunctions')
addpath('../ProductionFunctions')
addpath('../Economies')
addpath('../')

dir_out = '../../../data/generated/matlab/Calibration/';
dir.raw = '../../../data/raw/';
dir.generated = '../../../data/generated/';
%% prepare objects
% load calibrated parameters for skill distribution
par_cal = tools.csv2struct([dir_out, version_skill_dist,'_abilityDistThreeFixC_sol.csv']);
mom_cal = tools.csv2struct([dir_out, version_skill_dist,'_abilityDistThreeFixC_mom.csv']);

const = tools.csv2struct([dir.raw,'constants.csv']);
inc = tools.csv2struct([dir.generated,'r/mean_inc.csv']);

DensSkill_0 = abilityDistFixed(par_cal);
DensSkill_1 = abilityDistFixed(par_cal);

DensPart = ParticipationNormal(par_cal);

% utility function object
par_util.epsilon = const.epsilon;
Util = UtilityQuasiLinear(par_util);

% combine in SkillAndParticipation object
par_dist_0.dens_skill = DensSkill_0;
par_dist_0.dens_part = DensPart;
par_dist_0.util = Util;

par_dist_0.scale_inc = par_cal.scale_inc;

Dist_0 = SkillAndParticipation(par_dist_0);

par_dist_1 = par_dist_0;
par_dist_1.dens_skill = DensSkill_1;

Dist_1 = SkillAndParticipation(par_dist_1);

% Production function
par_pf.parse_inputs = false;

ProdFun = KrusellRobotsEach(par_pf);

% Economy object
par_econ.pf = ProdFun;
par_econ.util = Util;

par_econ.w_lb = 1e-6;
par_econ.w_ub = 1e4;

par_econ.lb = 1e-12 .* ones(1,6);
par_econ.ub = 1e6 .* ones(1,6);

par_econ.par_dist = par_cal;

par_econ_0 = par_econ;
par_econ_1 = par_econ;

% robot prices and taxes
% prices are not actually used, but not providing them leads to error
par_econ_0.q_B = 1;
par_econ_0.q_E = 1;
par_econ_0.q_S = 1;

par_econ_1.q_B = 0.4;
par_econ_1.q_E = 1;
par_econ_1.q_S = 1;

%% compute tax rates for capital

% tax rates on flows from Slavik & Yazici (2014) and net return on assets
% loaded from const

% net return is used to convert taxes on flow to taxes on stock
tax_equipment = const.tax_equipment_flow*const.net_return/(1+const.net_return);
tax_structures = const.tax_structures_flow*const.net_return/(1+const.net_return);

par_econ_0.dist = Dist_0;
par_econ_1.dist = Dist_1;

econ_0 = KrusellWithRobotsEachEconomy(par_econ_0);
econ_1 = KrusellWithRobotsEachEconomy(par_econ_1);

%% targets
targets.Y_M_0 = par_cal.Y_M_0;
targets.Y_R_0 = par_cal.Y_R_0;
targets.Y_C_0 = par_cal.Y_C_0;
targets.Y_M_1 = par_cal.Y_M_1;
targets.Y_R_1 = par_cal.Y_R_1;
targets.Y_C_1 = par_cal.Y_C_1;
targets.Y_B_0 = econ_0.q_B*(1+tax_equipment);
targets.Y_E_0 = econ_0.q_E*(1+tax_equipment);
targets.Y_S_0 = econ_0.q_S*(1+tax_structures);
targets.Y_E_1 = targets.Y_E_0; % only used if prices for capital fixed
targets.Y_S_1 = targets.Y_S_0; % only used if prices for capital fixed
targets.labor_income_share_0 = const.labor_share;
targets.M_share_in_labor_income_0 = inc.M_share;
targets.R_share_in_labor_income_0 = inc.R_share;
targets.labor_share_in_labor_plus_robot_income_0 = 0.9916; 
% above target not used, but keep to not mess things up
targets.K_B_rel_change = const.K_B_1_data/const.K_B_0_data - 1;
targets.alpha = const.KORV_alpha;
targets.rho = const.KORV_rho;
targets.sigma = const.KORV_sigma;
targets.Y_B_rel_change = -0.6; % works better with larger target
% (solution closer to true value)
targets.robot_share_capital = const.robot_share_capital;

varnames = {'L_M_0','L_R_0','L_C_0','K_B_M_0','K_B_R_0','K_B_C_0', ...
    'K_E_0','K_S_0',...
    'L_M_1','L_R_1','L_C_1','K_B_M_1','K_B_R_1','K_B_C_1', ...
    'kappa_M', ...
    'kappa_R', ...
    'kappa_C', ...
    'nu', ...
    'lambda', ...
    'mu', ...
    'eta_M', ...
    'eta_R', ...
    'eta_C', ...
    'upsilon', ...
    'A',...
    'alpha',...
    'rho',...
    'sigma',...
    'gamma',...
    'scale_L',...
    'tau_B',...
    'tau_E',...
    'tau_S'};


cal = CalibrateKrusellWithRobotsEachEconomy(econ_0, econ_1, targets, varnames,'num_nodes_f',100,'version',version_pf);
cal.fix_equipment_structures = true; % whether to keep equipment
% and structures fixed in
% period 1

smallNum = 1e-3;
smallNum_K_B = 1e-4;
bigNum = 1e4;

inp.K_B_M_0 = 86.5973;
inp.K_B_R_0 = 1.0000e-04;
inp.K_B_C_0 = 576.8583;
inp.K_E_0 = 1.2457e+03;
inp.K_S_0 = 762.8050;
inp.L_M_1 = 0.2433;
inp.L_R_1 = 1.6741;
inp.L_C_1 = 1.0223;
inp.K_B_M_1 = 986.8628;
inp.K_B_R_1 = 2.9071e-04;
inp.K_B_C_1 = 1.0404e+03;
inp.kappa_M = 3.2261;
inp.kappa_R = 2.8286;
inp.kappa_C = 1.1341;
inp.nu = 1.0055e+03;
inp.lambda = 0.6561;
inp.mu = 0.4568;
inp.eta_M = 0.7770;
inp.eta_R = 0.9990;
inp.eta_C = 0.7754;
inp.upsilon = 0.5249;
inp.A = 0.1196;
inp.alpha = 0.1142;
inp.rho = 0.5169;
inp.sigma = 1.4457;
inp.gamma = 1;
inp.scale_L = 255.8724; % scales labor in the prod. fct. (should
% include in table in paper?)
inp.tau_B = tax_equipment;
inp.tau_E = tax_equipment;
inp.tau_S = tax_structures;

% compute aggregate labor supplies
inp.L_M_0 = econ_0.L_M(par_cal.Y_M_0,par_cal.Y_R_0,par_cal.Y_C_0,par_cal);
inp.L_R_0 = econ_0.L_R(par_cal.Y_M_0,par_cal.Y_R_0,par_cal.Y_C_0,par_cal);
inp.L_C_0 = econ_0.L_C(par_cal.Y_M_0,par_cal.Y_R_0,par_cal.Y_C_0,par_cal);

inp.L_M_1 = econ_1.L_M(par_cal.Y_M_1,par_cal.Y_R_1,par_cal.Y_C_1,par_cal);
inp.L_R_1 = econ_1.L_R(par_cal.Y_M_1,par_cal.Y_R_1,par_cal.Y_C_1,par_cal);
inp.L_C_1 = econ_1.L_C(par_cal.Y_M_1,par_cal.Y_R_1,par_cal.Y_C_1,par_cal);

% lower bounds
lb.L_M_0 = inp.L_M_0;
lb.L_R_0 = inp.L_R_0;
lb.L_C_0 = inp.L_C_0;
lb.K_B_M_0 = smallNum_K_B;
lb.K_B_R_0 = smallNum_K_B;
lb.K_B_C_0 = smallNum_K_B;
lb.K_E_0 = smallNum;
lb.K_S_0 = smallNum;

lb.L_M_1 = inp.L_M_1;
lb.L_R_1 = inp.L_R_1;
lb.L_C_1 = inp.L_C_1;
lb.K_B_M_1 = smallNum_K_B;
lb.K_B_R_1 = smallNum_K_B;
lb.K_B_C_1 = smallNum_K_B;

lb.kappa_M = smallNum;
lb.kappa_R = smallNum;
lb.kappa_C = smallNum;
lb.nu = smallNum;
lb.lambda = smallNum;
lb.mu = smallNum;

lb.alpha = smallNum;
lb.rho = smallNum;
lb.sigma = smallNum;

lb.eta_M = smallNum;
lb.eta_R = smallNum;
lb.eta_C = smallNum;

lb.upsilon = smallNum;
lb.A = smallNum;
lb.gamma = 1;
lb.scale_L = 1e-3;

lb.tau_B = inp.tau_B;
lb.tau_E = inp.tau_E;
lb.tau_S = inp.tau_S;

% upper bounds
ub.L_M_0 = inp.L_M_0;
ub.L_R_0 = inp.L_R_0;
ub.L_C_0 = inp.L_C_0;
ub.K_B_M_0 = bigNum;
ub.K_B_R_0 = bigNum;
ub.K_B_C_0 = bigNum;
ub.K_E_0 = bigNum;
ub.K_S_0 = bigNum;

ub.L_M_1 = inp.L_M_1;
ub.L_R_1 = inp.L_R_1;
ub.L_C_1 = inp.L_C_1;
ub.K_B_M_1 = bigNum;
ub.K_B_R_1 = bigNum;
ub.K_B_C_1 = bigNum;
ub.kappa_M = bigNum;
ub.kappa_R = bigNum;
ub.kappa_C = bigNum;
ub.nu = bigNum;
ub.lambda = 1-smallNum;
ub.mu = 1-smallNum;

ub.alpha = 1-smallNum;
ub.rho = bigNum;
ub.sigma = bigNum;

ub.eta_M = 1-smallNum;
ub.eta_R = 1-smallNum;
ub.eta_C = 1-smallNum;

ub.upsilon = 1-smallNum;
ub.A = bigNum;
ub.gamma = 1;
ub.scale_L = 1e6;

ub.tau_B = inp.tau_B;
ub.tau_E = inp.tau_E;
ub.tau_S = inp.tau_S;
%% solve and save
sol = cal.calibrate(inp, lb, ub)
mom = cal.compute_moments(sol);

fprintf('nonlcon at solution\n')
[c,ceq] = cal.nonlcon_ssr(sol)


mom.revenue_requirement = mom_cal.income_tax_revenue + ...
    mom.tax_rev_equipment + mom.tax_rev_robots + mom.tax_rev_structures ...
    - par_cal.transfer * (1-mom_cal.participation_rate_0);

mom

tools.struct2csv(sol,[dir_out, version_pf,'_KrusellRobotsEach_sol.csv'])
tools.struct2csv(mom,[dir_out, version_pf,'_KrusellRobotsEach_mom.csv'])

end